Method for establishing the excitation force applied by the swell incident on a movable means of a wave energy system using a model of the drag force

ABSTRACT

The present invention is a method for real-time determination of the forces exerted by incident waves on a mobile part of a wave energy system. Models are constructed of the radiation force exerted on the mobile part and of the drag force exerted on the mobile part and a non-linear model of the wave energy system dynamics. The invention uses only measurements of the float kinematics (position, velocity and possibly acceleration) and of the force applied by a converter machine, which measurements are normally available on a wave energy system since they are used for control and supervision thereof. Determination of the excitation force exerted by incident waves on the mobile part uses the models, the measurements and an unscented Kalman filter.

CROSS-REFERENCE TO RELATED APPLICATIONS

Reference is made to International Application No. PCT/EP2019/050381,filed Jan. 9, 2019, which claims priority to French Patent ApplicationNo. 18/50.782, filed Jan. 31, 2018, the contents of which areincorporated herein by reference in their entirety.

BACKGROUND OF THE INVENTION Field of the Invention

The invention relates to conversion of wave energy to electrical energy.In particular, the invention concerns the determination of theexcitation force exerted by incident waves on a mobile part, notably forcontrolling a wave energy system.

Renewable energy resources have generated strong interest for someyears. They are clean, free and inexhaustible, which are major assets ina world facing the inexorable depletion of the available fossilresources and recognizing the need to preserve the planet. Among theseresources, wave energy which is a source relatively unknown amidst thosewidely publicized, such as wind or solar energy, contributes to thevital diversification of the exploitation of renewable energy sources.The devices, commonly referred to as “wave energy devices”, areparticularly interesting because they allow electricity to be producedfrom this renewable energy source (the potential and kinetic waveenergy) without greenhouse gas emissions. They are particularly wellsuited for providing electricity to isolated island sites.

Description of the Prior Art

For example, patent applications FR-2,876,751 and WO-2009/081,042describe devices intended to capture the energy produced by the seawater forces. These devices are made up of a floating support structureon which is a pendulum movably mounted with respect to the floatingsupport. The relative motion of the pendulum with respect to thefloating support is used to produce electrical energy by use of anenergy conversion machine (an electrical machine for example). Theconversion machine operates as a generator and as a motor. Indeed, inorder to provide torque or a force for driving the mobile part, power issupplied to the converter machine to bring it into resonance with thewaves (motor mode). On the other hand, to produce a torque or a forcethat withstands the motion of the mobile part, power is recovered viathe conversion machine (generator mode).

In order to improve the efficiency and therefore the cost-effectivenessof the devices for converting wave energy into electrical energy (waveenergy systems), it is interesting to estimate the excitation forceexerted by incident waves on the wave energy system.

The force applied by the waves to the devices for converting wave energyinto electrical energy or other forms of exploitable energy, commonlyreferred to as wave energy devices, cannot be directly measured undernormal operating conditions in an industrial context. For some waveenergy systems, it is possible to carry out specific tests which lockthe floating part of the machine and measure the force required tocounteract the action of the waves (i.e. to keep the float stationary)by use of one or more (force or torque) sensors arranged on theconversion machine, also referred to as power take-off (PTO) system,which for example are an electrical generator provided with a deviceallowing the oscillating motion of the mobile part to be transmitted.According to the common terminology, the quantity thus measured isreferred to as excitation force of the wave motion, by differentiatingit notably from the force (and the wave motion) generated by the verymotion of the float (radiation force). On the other hand, during normaloperation of the float, the same force sensors only measure the forceapplied by the PTO, and not by the excitation force of the incidentwaves.

In principle, it would be possible to calculate all the forces appliedby the waves by integrating pressure measurements provided by sensorsdistributed all over the surface, but this solution is expensive and notvery robust, it would therefore be difficult to envisage in anindustrial context.

There is little scientific work in the field of real-time estimation ofthe excitation force (or torque) of the waves. The document PeterKracht, Sebastian Perez-Becker, Jean-Baptiste Richard, and BorisFischer. “Performance Improvement of a Point Absorber Wave EnergyConverter by Application of an Observer-Based Control: Results from WaveTank Testing” is mentioned by way of example. IEEE Transactions onIndustry Applications, July 2015, describes a wave force estimationalgorithm based on a bank of independent harmonic oscillators and aLuenberger observer. The results obtained with this method show that theestimations are significantly delayed (phase shifted) in relation to thetrue excitation force, which is annoying notably for use within thecontext of wave energy system control. More generally, this approachconsiders an (a-priori selected) frequency range of fixed frequency forthe wave spectrum. For the method to work under realistic conditions,where the spectrum is time variant, a very large number of frequenciesneeds to be taken into account, so that this approach imposes a heavycomputational burden. The document Bradley A. Ling. “Real-TimeEstimation and Prediction of Wave Excitation Forces for Wave EnergyControl Applications”, published in ASME 2015, 34^(th) InternationalConference on Ocean, Offshore and Arctic Engineering, provides anapproach based on an extended Kalman filter for reconstructing thespectrum by considering the wave excitation force as a (single) sinewave whose parameters (amplitude, frequency and phase) vary over time.However, the method can only be effective for waves within a very narrowfrequency band.

Furthermore, patent application FR-2,973,448 corresponding to U.S. Pat.No. 9,261,070 describes a control method for oscillating pointconverters, One step estimates (the spectrum of) the wave force on thefloat (or of the torque on the mobile part of the converter), from a setof sine waves and a Luenberger observer. The method is similar to theone published in the aforementioned documents, and it involves a priorithe same drawbacks.

Another method is described in patent application FR-3,049,989(WO-2017/174,244). This method determines in real time the forcesexerted by incident waves on the mobile part, in order to adopt the bestadjustments for electrical energy production strategies. The method isbased on the construction of a model of the radiation force exerted onthe mobile part and of a model of the wave energy system dynamics. Thismethod only uses measurements of the float kinematics (position,velocity and possibly acceleration) and of the force applied by theconverter machine, which measurements are normally available on a waveenergy system since they are used for control and supervision thereof.Thus, the used models allow estimating the force exerted by the wavesfor any wave frequency range, while keeping a computation time suited toreal-time determination and control.

This method performs satisfactorily for “linear” wave energy systems,that is whose dynamic behavior is described in a sufficiently precisemanner by linear differential equations. However, this method does notallow the wave force to be determined precisely for “non-linear” waveenergy systems (with non-linear differential equations), notably flaptype wave energy systems or wave energy systems of flat submergedpressure differential converter type, or wave energy systems of surfacefloating rotating mass device type.

SUMMARY OF THE INVENTION

In order to overcome these drawbacks, the present invention concerns amethod for real-time determination of the forces exerted by incidentwaves on a mobile part of a wave energy system, in order to adopt thebest adjustments for electrical energy recovery strategies. The methodis based on the construction of a model of the radiation force exertedon the mobile part, a model of the drag force exerted on the mobile partand of a non-linear model of the wave energy system dynamics. Theinvention uses only measurements of the float kinematics (position,velocity and possibly acceleration) and of the force applied by theconversion machine, which measurements are normally available on a waveenergy system since they are used for control and supervision thereof.Determination of the excitation force exerted by incident waves on themobile part uses these models, the measurements and an unscented Kalmanfilter. It notably integrates a model of the viscous drag force and itallows, by use of the unscented Kalman filter, estimation of the forceexerted by incident waves in a precise manner for all wave energysystems, in particular the “non-linear” wave energy systems. The methodaccording to the invention is suited for any wave frequency range, whileproviding a a computation time suited for real-time determination andcontrol.

The present invention relates to a method of determining an excitationforce exerted by the waves on a mobile part of a wave energy system. Thewave energy system converts the wave energy into to electrical energythrough the mobile part cooperating with a conversion machine, themobile part moving with respect to the conversion machine under theaction of the waves. For this method, the following steps are carriedout:

-   -   a) measuring the position and the velocity of the mobile part;    -   b) measuring the force u exerted by the converter machine on the        mobile part;    -   c) constructing a model of the radiation force τ_(rad) exerted        on the mobile part, the radiation force model relating the        radiation force τ_(rad) to the velocity of the mobile part;    -   d) constructing a model of the drag force τ_(d) exerted on the        mobile part, the drag force model relating the drag force τ_(d)        to the velocity of the mobile part;    -   e) constructing a dynamic model of the wave energy system        relating the excitation force τ_(w) exerted by the incident        waves on the mobile part to the position of the mobile part, to        the velocity of the mobile part, to the force u exerted by        conversion machine on the mobile part, to the radiation force        τ_(rad) exerted on the mobile part and to the drag force τ_(d)        exerted on the mobile part; and    -   f) determining the excitation force τ_(w) exerted by incident        waves on the mobile part using the dynamic model, the radiation        force model, the drag force τ_(d) model, the measured position        and velocity, and the measured force u exerted by the conversion        machine on the mobile part, and using an unscented Kalman filter        constructed from a random walk model of the excitation force        exerted by the incident waves on the mobile part.

According to an embodiment, the dynamic model of the wave energy systemis constructed by use of an equation:

(t)=τ_(hd)(t)+τ_(rad)(t)+τ_(d)(t)+τ_(w)(t)−u(t). with

being the total moment of inertia of the mobile part [2]. δ(t) being theangle of rotation of the mobile part with respect to the equilibriumposition, with

(t) being the angular acceleration of the mobile part and

(t) the angular velocity of the mobile part, τ_(hd)(t) being thehydrostatic restoring moment. τ_(rad)(t) being the radiation moment,τ_(d)(t) being the drag moment. τ_(w)(t) being the wave excitationmoment and u(t) being the moment exerted by the conversion machine onthe mobile part.

According to an implementation, hydrostatic restoring moment τ_(hd)(t)is determined by use of a formula:

τ_(hd)(t) = −K δ(t)

where K is the hydrostatic stiffness coefficient.

According to an aspect, the radiation force model is constructed by useof an equation of the type;

τ_(rad)(t)=−

(t)−τ_(r)(t), with

being the added moment of inertia at infinite high frequency; andτ_(r)(t)=

h(t−s)

(s)ds=h(t)*

(t), with h being the impulse response relating the velocity of themobile part to the radiation damping and

(t) being the angular velocity of the mobile part.

Advantageously, the drag force model is constructed by use of anequation of the type.

${\tau_{d}(t)} = {\beta \; {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$

where β is the drag coefficient and

(t) is the angular velocity of the mobile part.

Alternatively, the drag force model F_(d) can be written by use of anequation of the type:

${F_{d}(t)} = {\beta \; {\overset{.}{z}(t)}{{z(t)}}}$

with β being the drag coefficient and

being the velocity of the mobile part.

According to an implementation. the excitation force τ_(w) exerted bythe incident waves on the mobile part is determined by carrying out thefollowing steps:

-   -   i) initializing k=0, state vector {circumflex over        (x)}_(α)(0|0)=m(0) and the state of the covariance matrix        P(0|0)=    -   ii) at any time k, acquiring the position and velocity        measurements of the mobile part γ(k)=[δ(k)        (k)] and the measurement of the force exerted by the converter        machine on the mobile part u(k); and    -   iii) at any time k, determining the excitation force exerted by        the incident waves on the mobile part        _(w)(k) by use of the following equations:

     τ̂_(w)(k) = [0  1]x̂_(a)(kk)     P_(x)(kk) = P_(x)(kk − 1) − KP_(y)(kk − 1)K^(T)     K = P_(xy)(kk − 1)P_(y)(kk − 1)⁻¹     x̂(k) = x̂_(a)(kk − 1) + K(ŷ(k) − ŷ(kk − 1))${P_{xy}\left( {k{k - 1}} \right)} = {\sum\limits_{j = 0}^{2n}\; {{W_{j}^{e}\left( {{{\hat{x}}_{j}\left( {k{k - 1}} \right)} - {x_{a}\left( {k{k - 1}} \right)}} \right)}\left( {{{\hat{y}}_{j}\left( {k{k - 1}} \right)} - {y\left( {k{k - 1}} \right)}} \right)^{T}}}$$\mspace{76mu} {{{\hat{y}\left( {k{k - 1}} \right)} = {\sum\limits_{j = 0}^{2n}\; {W_{j}^{m}{{\hat{y}}_{j}\left( {k{k - 1}} \right)}}}},{{P_{y}\left( {k{k - 1}} \right)} = {{\sum\limits_{j = 0}^{2n}\; {{W_{j}^{e}\left( {{{\hat{y}}_{j}\left( {k{k - 1}} \right)} - {y\left( {k{k - 1}} \right)}} \right)}\left( {{{\hat{y}}_{j}\left( {k{k - 1}} \right)} - {y\left( {k{k - 1}} \right)}} \right)^{T}}} + R}}}$     ŷ_(j)(kk − 1) = C_(a)x̂_(j)(kk − 1)$\mspace{76mu} {{{\hat{x}\left( {k{k - 1}} \right)} = {\sum\limits_{j = 0}^{2n}\; {W_{j}^{m}{{\hat{x}}_{j}\left( {k{k - 1}} \right)}}}},{{P_{x}\left( {k{k - 1}} \right)} = {{\sum\limits_{j = 0}^{2n}\; {{W_{j}^{e}\left( {{{\hat{x}}_{j}\left( {k{k - 1}} \right)} - {x_{a}\left( {k{k - 1}} \right)}} \right)}\left( {{{\hat{x}}_{j}\left( {k{k - 1}} \right)} - {x_{a}\left( {k{k - 1}} \right)}} \right)^{T}}} + Q}}}$x̂_(j)(kk − 1) = A_(a)x_(j)(k − 1) + f_(a)(x_(j)(k − 1)) + B_(a)u(k − 1), j = 0, 1, … , 2n$\mspace{76mu} {{{x_{0}\left( {k - 1} \right)} = {{\hat{x}}_{a}\left( {{k - 1}{k - 1}} \right)}},\mspace{76mu} {{x_{i}\left( {k - 1} \right)} = {{{\hat{x}}_{a}\left( {{k - 1}{k - 1}} \right)} + {\sqrt{n + \lambda}{S_{i}\left( {k - 1} \right)}}}},{i = 1},2,\ldots \;,n}$$\mspace{76mu} {{{x_{i + n}\left( {k - 1} \right)} = {{{\hat{x}}_{a}\left( {{k - 1}{k - 1}} \right)} - {\sqrt{n + \lambda}{S_{i}\left( {k - 1} \right)}}}},{i = 1},2,\ldots \;,n}$

where S_(i)(k−1) is the i-th column of the matrix square root ofP_(x)(k−1|k−1), that is

P_(x)(k|k − 1) = S(k − 1)^(T)S(k − 1).

x_(α)(k) being the unknown state vector

${{x_{c}(k)} = \begin{bmatrix}{x(k)} \\{\tau_{w}(k)}\end{bmatrix}},$

matrices of the state-space realization. P being the covariance matrixof the state vector, Q and R calibration matrices.

According to an embodiment, the wave energy system is controlleddepending on the determined excitation force τ_(w) exerted by theincident waves on the mobile part.

Advantageously, the wave energy system is a wave energy system of a flatsubmerged pressure differential converter construction or of surfacefloating rotating mass device type, or a flap type wave energy system.

Furthermore, the invention relates to a wave energy system. It comprisescontrol, notably computer control, implementing a method according toone of the above features.

BRIEF DESCRIPTION OF THE DRAWINGS

Other features and advantages of the method according to the inventionwill be clear from reading the description hereafter, with reference tothe accompanying figures wherein:

FIG. 1 illustrates a flap construction wave energy system according toan embodiment of the invention;

FIG. 2 illustrates an estimation of the force exerted by the waves byuse of a method according to the prior art for a regular wave;

FIG. 3 illustrates an estimation of the force exerted by the waves byuse of a method according to an embodiment of the invention for aregular wave; and

FIG. 4 illustrates an estimation of the force exerted by the waves byuse of a method according to an embodiment of the invention for anirregular wave.

DETAILED DESCRIPTION OF THE INVENTION

The present invention relates to a method of determining the excitationforce exerted by incident waves on a mobile part of a wave energysystem, which is also referred to as wave excitation force. A waveenergy system is a system that converts wave energy to recoverableenergy, in particular electrical energy. A wave energy system generallycomprises a mobile part, also referred to as flap, pendulum or float,which has an oscillating motion under the action of the waves. Themobile part cooperates with a conversion machine, also referred to as apower take-off (PTO) system, which comprises in most cases an electricalgenerator coupled to a device allowing adaption of the transmission ofthe oscillating motion, in order to convert the motion of the mobilepart into recoverable energy. In some cases, the conversion machine canact as a motor by generating a force on the mobile part. Indeed, inorder to recover power via the conversion machine, a torque or a forceis produced opposing the motion of the mobile (generator mode).Furthermore, if the conversion machine allows, power can be supplied tothe conversion machine to provide torque or a force for driving themobile part in order to bring it into resonance with the waves (motormode).

The method according to the invention is suited for any type of waveenergy system with at least one mobile part, as described for example inpatent application FR-2,973,448 corresponding to U.S. Pat. No.9,261,070. The control method according to the invention can also beapplied to a wave energy system belonging to the category of wave energysystems with oscillating water columns (OWC). However, the methodaccording to the invention is particularly suited for a non-linear waveenergy system, for example of flap construction, of flat submergedpressure differential converter construction or of surface floatingrotating mass construction.

FIG. 1 schematically illustrates a non-limitative example of a flapdesign wave energy system 1. Wave energy system 1 comprises a mobilepart 2 in a form of a substantially vertical flap (the wave motion isschematically represented by two curved lines and the direction of thewaves is shown by arrow H). Flap 2 has a submerged part and an emergedpart. Mobile part 2 cooperates, through a hinge 4, with a convertermachine based on an electrical generator 3, which is in this case arotary generator. Hinge 4 allows a rotary oscillating motion of flap 2about a horizontal axis with respect to a support 5 that is stationarywith respect to the seabed 6.

In the rest of the description, the terms waves, sea waves and wavemotion are considered to be equivalent.

Furthermore, in the description, the term force designates a stress or atorque. Similarly, the terms position, velocity and accelerationdesignate both “linear” and “angular” values. The “linear” values can beassociated with the stress and the “angular” values can be associatedwith the torque. In the rest of the description, only the case relativeto torques is illustrated, but the stress-related case can be deduced bytransposing the equations to an orthogonal reference frame.

Moreover, for a better understanding, the various models are representedin one dimension. The method according to the invention is suited formodels in several dimensions, for systems whose motion has severaldegrees of freedom.

The method according to the invention comprises the following steps:

-   -   1. Measurement of the mobile part's position and velocity    -   2. Measurement of the force exerted by the converter machine (or        PTO)    -   3. Construction of the radiation force model    -   4. Construction of the drag force model    -   5. Construction of the dynamic model    -   6. Determination of the incident wave excitation force    -   7. (optional step) Control of the wave energy system.    -   Steps 3 and 4 can be carried out in this order, in the reverse        order or simultaneously.

The steps of the method can be carried out by a computer equipment (orcalculator). This equipment can comprise data processing and,advantageously, data storage.

-   -   The data processors are configured to implement:        -   measurement of the mobile part's position and velocity,        -   measurement of the force exerted by the converter machine,        -   construction of the radiation force, drag force and dynamic            models, and        -   determination of the incident wave excitation force.

1—Measurement of the Mobile Part's Position and Velocity

This step measures the position and the velocity of the mobile part. Theposition corresponds to the motion (distance or angle for example) withrespect to the position of equilibrium of the mobile part. Thesemeasurements can be performed using sensors, generally present on a waveenergy system for at least one of control and supervision thereof.

According to an implementation of the invention, this step can alsoinclude measuring the acceleration of the mobile part, which can be usedto estimate the velocity, or directly in the models used by the methodaccording to the invention. For example, acceleration can be measuredusing an accelerometer arranged on the mobile part.

2—Measurement of the Force Exerted by the Converter Machine (PTO)

The force (the stress or possibly the torque) exerted by the PTOconverter machine on the mobile part is measured in this step. Thismeasurement can be performed using a sensor, which can be a force sensoror a torque sensor. This type of sensor is often installed or it can bereadily installed in wave energy systems, for at least one of controland supervision thereof. Alternatively, the measurement can be replacedby an estimation performed from the force (or torque) setpoint sent tothe PTO.

For the wave energy system example illustrated in FIG. 1, a torquesensor can be arranged at hinge 4 (or electrical generator 3).

3—Construction of the Radiation Force Model

This step constructs a model of the radiation force exerted on themobile part. According to the linear wave theory (as described forexample in the document Falnes J, Kurniawan A. “Fundamental formulae forWave-Energy Conversion”. R. Soc. open sci. 2: 140305, 2005,http://dx.doi.org/10.1098/rsos.140305), the radiation force results fromthe oscillation of a submerged body (therefore is depends on the motionof the mobile part), while the excitation force, resulting from the verypresence of a body in the water, does not depend on the motion of thesubmerged body, but is dependent on the incident wave. In the absence ofincident wave, the radiation force damps the residual oscillation of thesubmerged body and eventually stops it. It is important to note that,although the linear theory allows relating the excitation force to theelevation of the incident wave through a linear model (in the frequencyor time domain), in practice it cannot be used to calculate theexcitation force online, even though it is possible to measure theelevation of the wave at the center of gravity of the mobile part asrequired by the theory. Indeed, the linear relation between waveelevation and excitation force is non-causal, which means that theexcitation force at a given time cannot be calculated without knowingthe wave elevation in future times (on the other hand, calculation canbe carried out offline, once the wave has passed). In a real-timecontrol context, the excitation force can therefore only be consideredas a totally unknown exogenous force acting on the float. On the otherhand, still according to the linear wave theory, the radiation force isrelated to the motion of the float, and more precisely to theacceleration and the velocity thereof, by a causal linear model (in thefrequency or time domain). It can therefore be calculated online usingthe current acceleration measurements and the current and past velocitymeasurements.

According to an implementation of the invention, the radiation forcemodel

According to an implementation of the invention, the radiation forcemodel τ_(rad)(t) is constructed by use of an equation:

${\tau_{rad}(t)} = {{{- I_{\infty}}{\overset{¨}{\delta}(t)}} - {\tau_{r}(t)}}$

with:τ_(r)(t)=

h(t−s)

(s)ds=h(t)

(t) the component of the radiation force τ_(rad) that depends on the(current and past) velocity of the mobile part, which can be referred toas radiation damping:

is the angular acceleration of the mobile part:

is the added moment of inertia at infinite high frequency, which can beobtained using BEM (Boundary Element Method) calculation codes, such asWAMIT®[WAMIT, USA] or Nemah® (Ecole Centrale de Nantes, France), fromthe geometry of the mobile part:

is the angular velocity of the mobile part, andh is the impulse response relating the velocity of the mobile part tothe radiation damping, whose coefficients are obtained from hydrodynamicparameters of the mobile part calculated with the same BEM calculationcodes.

The construction of this model allows determination of radiation forceof any time, with a limited computation time. Thus, determination of theforce exerted by the waves can be determined at any time with a shortcomputation time.

4—Construction of the Drag Force Model

This step constructs a model of the drag force exerted on the mobilepart. The drag force corresponds to the viscous friction on the mobilepart of the wave energy system. It is a non-linear model that can dependon the velocity of the mobile part. The drag force (due to viscousfriction) is often considered to be negligible for wave energy systemsof oscillating point converter type and it is generally excluded fromthe modeling thereof. This is not the case, however, for flapconstruction wave energy systems or for machines such as flat submergedpressure differential converters or surface floating rotating massdevices, among others.

According to an embodiment of the invention, for a wave energy systemhaving a rotational motion, the drag force model τ_(d) can be writtenwith an equation:

${\tau_{d}(t)} = {\beta \; {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$

with β being the drag coefficient, this coefficient can be determined byexperimental tests of the wave energy system or CFD type (ComputationalFluid Dynamics) numerical simulations; and

being the angular velocity of the mobile part.

According to a variant embodiment of the invention, for a wave energysystem having a translational motion, the drag force model F_(d) can bewritten with an equation:

${F_{d}(t)} = {\beta \; {\overset{.}{z}(t)}{{\overset{.}{z}(t)}}}$

with β being the drag coefficient which can be determined byexperimental tests of the wave energy system or CFD type (ComputationalFluid Dynamics) numerical simulations, and

being the velocity of the mobile part.

The construction of this model enables precise determination of the dragforce. Thus, determination of the force exerted by the waves on themobile part can be precise because it does not disregard the drag force.

5—Construction of the Dynamic Model of the Wave Energy System

This step constructs a dynamic model of the wave energy system. Adynamic model is understood to be a model that relates the excitationforce exerted by the incident waves on the mobile part, the radiationforce exerted on the mobile part, the hydrostatic restoring forceexerted on the mobile part, the force exerted by the converter machineon the mobile part, to the position and the velocity of the mobile part.This type of model allows obtaining results representative of thebehavior of the wave energy system if the motions are not too large.

Advantageously, the dynamic model is obtained by applying thefundamental principle of dynamics (Newton's second law) to the mobilepart.

According to an embodiment of the invention wherein the stresses areconsidered, the dynamic model of the wave energy system can beconstructed with an equation:

${I_{eq}{\overset{¨}{\delta}(t)}} = {{\tau_{hd}(t)} + {\tau_{{ra}\; d}(t)} + {\tau_{d}(t)} + {\tau_{w}(t)} - {u(t)}}$

with:

being the moment of inertia of the mobile part,

being the angular acceleration of the mobile part,τ_(w) being the excitation torque exerted by the incident waves on themobile part,τ_(rad) being the radiation torque exerted on the mobile part,τ_(hd) being the hydrostatic restoring torque exerted on the mobilepart,τ_(d) being the drag torque exerted an the mobile part, andu being the torque exerted by the converter machine on the mobile part.

This model conveys a rotational motion about a horizontal axis (typicalfor the wave energy system of FIG. 1). This model is derived from thelinear wave theory.

According to a first variant embodiment, the hydrostatic restoring forceexerted on the mobile part can be approximated by a linear function ofposition z defined with respect to the equilibrium position. In thiscase, the hydrostatic restoring force can be written by a function ofthe type: τ_(hd)(t)=−Kδ(t), with δ being the angular position of themobile part defined with respect to its equilibrium position and K thehydrostatic stiffness coefficient. Thus, the hydrostatic restoring forcecan be calculated from a simple model if the measurement of position δis available. This function is particularly well-suited for smalldisplacements δ.

According to an implementation of the invention wherein the stresses areconsidered, the dynamic model of the wave energy system can beconstructed with an equation:

${M{\overset{¨}{z}(t)}} = {{F_{ex}(t)} + {F_{hd}(t)} + {F_{r\; {ad}}(t)} + {F_{d}(t)} - {F_{u}(t)}}$

with:M being the mass of the mobile part,

being the acceleration of the mobile part,

being the excitation force exerted by the incident waves on the mobilepart,F_(rad) being the radiation force exerted on the mobile part,

being the hydrostatic restoring force exerted on the mobile part,F_(d) being the drag force exerted on the mobile part, andF_(u) being the force exerted by the converter machine on the mobilepart.

This model conveys a vertical translational motion (typical of floatshaving a heave motion). This model is derived from the linear wavetheory.

This model is the “mirror” of the model when the torques are considered;the various terms of the model are similar in nature.

According to an embodiment of the invention, at least one of thehydrostatic restoring force and the hydrostatic restoring torque can beapproximated by a linear function or by a piecewise affine function.

6—Determination of the Excitation Force Exerted by the Incident Waves

In this step, real-time determination of the excitation force exerted bythe incident waves on the mobile part is performed by means of:

-   -   the mobile part's position and velocity (and possibly        acceleration measurements determined in step 1;    -   the measurement of the force exerted by the converter machine        PTO on the mobile part determined in step 2    -   the radiation force model determined in step 3;    -   the drag force model determined in step 4; and    -   the wave energy system dynamic model determined in step 5.

According to the invention, the excitation force exerted by the incidentwaves on the mobile part is determined using an observer based on anunscented Kalman filter (UKF) constructed from a random walk model ofthe excitation force exerted by the incident waves on the mobile part.The unscented Kalman filter allows to take accounting for thenon-linearities of the models, in particular the drag force model.

It is noted that a state observer, or state estimator, is, in automationand systems theory, an extension of a model represented as a state-spacerepresentation. When the state of the system is not measurable, anobserver allowing the state to be reconstructed from a model isconstructed.

A UKF filter is based on the “unscented” transformation theory, whichallows an estimator to be obtained for a non-linear system withoutrequiring prior linearization for application to the filter. The UKFfilter uses a statistical state distribution that is propagated throughthe non-linear equations. Such a filter has the advantage of providingestimation stability and therefore robustness.

For this embodiment, what is known prior to this step thus is:

-   -   the mobile part's positon δ and velocity        measurements,    -   the measurement of force u exerted by the conversion machine PTO        on the mobile part,    -   the radiation force model τ_(rad)(t)=−        (t)−τ_(r)(t), with

${{\tau_{r}(t)} = {{\int_{0}^{t}{{h\left( {t - z} \right)}{\overset{.}{\delta}(s)}{ds}}} = {{h(t)}*{\overset{.}{\delta}(t)}}}},$

-   -   the drag force model τ_(d)(t)=β        (t)|        (t)|, and    -   the wave energy system dynamic model

${I_{eq}{\overset{¨}{\delta}(t)}} = {{\tau_{hd}(t)} + {\tau_{{ra}\; d}(t)} + {\tau_{d}(t)} + {\tau_{w}(t)} - {{u(t)}.}}$

In this approach, the problem of estimation of the wave excitation forceis transformed into a conventional state estimation problem (that can besolved with an unscented Kalman filter), by expressing the dynamics ofthe wave excitation force by a random walk model. The main advantage ofthis method is the consideration of uncertainties allowing accountingfor of the measurement noises and the unmodelled dynamics.

By replacing in the equation that describes the motion of the mobilepart

${I_{eq}{\overset{¨}{\delta}(t)}} = {{\tau_{hd}(t)} + {\tau_{r\; {ad}}(t)} + {\tau_{d}(t)} + {\tau_{w}(t)} - {u(t)}}$

the expressions for the hydrostatic restoring force (with the linearmodel), the radiation force model and the drag force model

τ_(hd)(t) = −K δ(t)${\tau_{\; {{ra}\; d}}(t)} = {{{- I_{\infty}}{\overset{¨}{\delta}(t)}} - {\tau_{r}(t)}}$${\tau_{d}(t)} = {\beta \; {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$

the following non-linear model is obtained:

${{\left( {I_{eq} + I_{\infty}} \right){\overset{¨}{\delta}(t)}} + {\beta {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}} + {K\; {\delta (t)}}} = {{\tau_{w}(t)} + {\tau_{r}(t)} - {{u(t)}.}}$

This equation can be put in state form, by defining

$\left\{ {\begin{matrix}{{x_{1}(t)} = {\delta (t)}} \\{{x_{2}(t)} = {\delta (t)}}\end{matrix}\quad} \right.$

which allows the previous model to be written in form of state-spacerepresentation

$\left\{ {\begin{matrix}{{{{\hat{x}}_{1}(t)} = {x_{2}(t)}}\mspace{526mu}} \\{{{\hat{x}}_{1}(t)} = {\frac{1}{I_{eq} + I_{oc}}\left( {{- {{Kx}_{1}(t)}} - {\beta \; {x_{2}(t)}{{x_{2}(t)}}} - {\tau_{r}(t)} + {\tau_{w}(t)} - {u(t)}} \right)}}\end{matrix}\quad} \right.$

This system of equations contains the integral term

${\tau_{r}(t)} = {{\int_{0}^{t}{{h\left( {t - s} \right)}{\overset{.}{\delta}(s)}{ds}}} = {{h(t)}*{\overset{.}{\delta}(t)}}}$

that can be considered to be a linear system, with

(t) as the input and τ_(r)(t) as the output. With Prony's method, thissystem can then be transformed into the equivalent state-spacerepresentation:

$\left\{ {\begin{matrix}{{{\hat{x}}_{r}(t)} = {{A_{r}{x_{r}(t)}} + {B_{r}{\delta (t)}}}} \\{{\tau_{r}(t)} = {{C_{r}{x_{r}(t)}} + {D_{r}{\delta (t)}}}}\end{matrix}{or}\left\{ \begin{matrix}{{{\hat{x}}_{r}(t)} = {{A_{r}{x_{r}(t)}} + {B_{r}{x_{2}(t)}}}} \\{{\tau_{r}(t)} = {{C_{r}{x_{r}(t)}} + {D_{r}{x_{2}(t)}}}}\end{matrix} \right.} \right.$

where x_(r) is an internal state [inaccessible] with no particularphysical meaning and (

) are the matrices of the state-space realization.

By combining the two state-space representations, it is obtained:

$\left\{ {\begin{matrix}{{{{\overset{.}{x}}_{1}(t)} = {x_{2}(t)}}\mspace{655mu}} \\{{{\overset{.}{x}}_{1}(t)} = {\frac{1}{I_{eq} + I_{oc}}\left( {{- {{Kx}_{1}(t)}} - {D_{r}{x_{2}(t)}} - {C_{r}{x_{r}(t)}} - {\beta \; {x_{2}(t)}{{x_{2}(t)}}} + {\tau_{w}(t)} - {u(t)}} \right)}} \\{{{{\overset{.}{x}}_{r}(t)} = {{B_{r}{x_{2}(t)}} + {A_{r}{x_{r}(t)}}}}\mspace{529mu}}\end{matrix}\quad} \right.$

or, in an equivalent manner

$\left\{ {\begin{matrix}{{\overset{.}{x}(t)} = {{A_{c}{x(t)}} + {f_{c}\left( {x(t)} \right)} + {B_{c}\left( {{\tau_{w}(t)} - {u(t)}} \right)}}} \\{{{y(t)} = {C_{c}{x(t)}}}\mspace{290mu}}\end{matrix}\quad} \right.$

where

x(t) = [x₁(t)  x₂(t)  x_(r)^(T)(t)]^(T) $y = \begin{bmatrix}{x_{1}(t)} \\{x_{2}(t)}\end{bmatrix}$

and

${A_{c} = \begin{bmatrix}0 & 1 & 0 \\{- \frac{K}{I_{eq} + I_{oc}}} & {- \frac{D_{r}}{I_{eq} + I_{oc}}} & {- \frac{C_{r}}{I_{eq} + I_{oc}}} \\0 & B_{r} & A_{r}\end{bmatrix}},{{f_{c}\left( {x(t)} \right)} = \begin{bmatrix}0 \\{{- \frac{\beta}{I_{eq} + I_{oc}}}{x_{2}(t)}{{x_{2}(t)}}} \\0\end{bmatrix}}$ $\mspace{76mu} {{B_{c} = \begin{bmatrix}0 \\\frac{1}{I_{eq} + I_{oc}} \\0\end{bmatrix}},{C_{c} = \begin{bmatrix}1 & 0 & 0 \\0 & 1 & 0\end{bmatrix}}}$

Because of the term ƒ_(c)(x(t)), the system is non-linear. The systemhas two inputs, τ_(w)(t) and u(t), and two outputs, x₁(t)=δ(t)=δ(t) andx₂(t)=

(t). Input τ_(w)(t) is not measurable and it is unknown. The problem tobe solved is to estimate it from the measured quantities u(t) [momentapplied by the converter machine PTO], δ(t) and

(t) [angular position and velocity of the mobile part].

If the principal motion of the wave energy system were a translationalmotion and the problem was to estimate the wave excitation force, thesame developments would apply by replacing angular position and velocitywith positions and velocities, and all the moments with forces.

To carry out the estimation, the above state system is first discretized(because the measurements are sampled and the estimation algorithm isexecuted by a calculator) using the Euler method, which yields, for agiven sampling period

:

$\left\{ {\begin{matrix}{{x\left( {k + 1} \right)} = {{A_{d}{x(k)}} + {f_{d}\left( {x(k)} \right)} + {B_{d}{\tau_{w}(k)}} - {B_{d}{u(k)}}}} \\{{{y(k)} = {C_{d}{x(k)}}}}\end{matrix}\quad} \right.$

where A_(d)=

ƒ_(d)(x(k))=T_(s)ƒ_(c)(x(k)), B_(d)=

C_(d)=

and

is the identity matrix of appropriate dimensions.

To estimate excitation moment τ_(w)(k), it is considered as a state, byintroducing a mathematical model that relates τ_(w)(k) and τ_(w)(k+1),in this case:

τ_(w)(k + 1) = τ_(w)(k) + ɛ_(m)(k)

where

(k) describes the variation of τ_(w)(k) and is considered to be a randomnumber. In other words, this model implies that, of any time k, theexcitation moment strays by a random step (quantity) from its previousvalue, and that these steps are independently and identically sizedistributed.

This random walk model of the excitation force is coupled with a morerealistic model of the mobile part dynamics:

$\left\{ {\begin{matrix}{{x\left( {k + 1} \right)} = {{A_{d}{x(k)}} + {f_{d}\left( {x(k)} \right)} + {B_{d}{\tau_{w}(k)}} - {B_{d}{u(k)}} + {\epsilon_{x}(k)}}} \\{{{y(k)} = {{C_{d}{x(k)}} + {v(k)}}}\mspace{374mu}}\end{matrix}\quad} \right.$

where

(k) represents the unmodelled dynamics (PTO friction, hydrostaticnon-linearity, etc.) and v(k) describes the measurement noise corruptingthe float position and velocity measurements.

By combining the random walk model of the excitation force and thenon-linear dynamics of the mobile part, the augmented system isobtained:

$\left\{ {\begin{matrix}{{x\left( {k + 1} \right)} = {{A_{d}{x(k)}} + {f_{d}\left( {x(k)} \right)} + {B_{d}{\tau_{w}(k)}} - {B_{d}{u(k)}} + {\epsilon_{x}(k)}}} \\{{{w\left( {k + 1} \right)} = {{\tau_{w}(k)} + {\epsilon_{m}(t)}}}\mspace{335mu}} \\{{{y(k)} = {{C_{d}{x(t)}} + {v(k)}}}\mspace{380mu}}\end{matrix}\quad} \right.$

or, in an equivalent manner:

$\left\{ {\begin{matrix}{{x_{a}\left( {k + 1} \right)} = {{A_{a}{x_{a}(k)}} + {f_{a}\left( {x_{a}(k)} \right)} + {B_{a}{u(k)}} + {\epsilon (k)}}} \\{{{y(k)} = {{C_{a}{x_{a}(k)}} + {v(k)}}}}\end{matrix}\quad} \right.$

where

${{x_{a}(k)} = \begin{bmatrix}{x(k)} \\{\tau_{w}(k)}\end{bmatrix}},{{f_{a}\left( {x_{a}(k)} \right)} = \begin{bmatrix}{f_{d}\left( {x(k)} \right)} \\0\end{bmatrix}},{{\epsilon (k)} = \begin{bmatrix}{ɛ_{x}(k)} \\{ɛ_{m}(k)}\end{bmatrix}}$

and

${A_{a} = \begin{bmatrix}A_{d} & B_{d} \\0 & 1\end{bmatrix}},{B_{a} = \begin{bmatrix}{- B_{d}} \\0\end{bmatrix}},{C_{a} = \left\lbrack {C_{d}\mspace{14mu} 0} \right\rbrack}$

Thus, the problem to be estimated, τ_(w)(k), becomes a state estimationproblem.

One way of estimating the unknown state vector x_(α)(k), which is byconsidering information on

(k) and v(k), applies the Kalman filter

algorithm. In this proposed method, the unscented Kalman filter (UKF) isused for handling the non-linearity of the system. The UKF is generallymore robust and precise than the extended Kalman filter

, which deals with the non-linearity by linearizing it. Furthermore, inthis case, the presence of the drag term, whose derivative is notcontinuous, makes the EKF inapplicable.

Like the EKF, the UKF performs the estimation in two steps, stateprediction and measurement correction, except that these two steps arepreceded by a prior step for “sigma points” calculation. The sigmapoints are a set of samples calculated so as to be able to exactlypropagate the mean and variance information in the space of a non-linearfunction.

According to an implementation of the invention, the followinghypotheses can be adopted:

-   -   the initial state x_(α)(0) is a random vector of mean        m(0)=E[x_(α)(0)] and of covariance        P(0)=E[(x_(α)(0)−m(0))(x_(α)(0)−m(0))^(T)].    -   (k) and v(k) are Gaussian noises with covariance matrices Q and        R respectively.        as well as the following notations:    -   {circumflex over (x)}_(α)(k|k−1) is the estimation of x_(α)(k)        from measurements up to the time k−1, i.e. γ(k−1), γ(k−2), . . .        and u(k−1), u(k−2), . . .    -   {circumflex over (x)}_(α)(k|k) is the estimation of x_(α)(k)        from measurements up to the time k, i.e. γ(k), γ(k−1) . . . and        u(k), u(k−1), . . .    -   (k|k−1) is the covariance matrix of x_(α)(k) from measurements        up to the time k−1, i.e. γ(k−1), γ(k−2), . . . and u(k−1),        u(k−2), . . .    -   (k|k) is the covariance matrix of x(k) from measurements up to        the time k. i.e. γ(k), γ(k−1) . . . and u(k), u(k−1), . . .        For this implementation, the three steps in the UKF method can        be:        1. Calculation of the sigma points

Let:

${W_{0}^{m} = \frac{\lambda}{n + \lambda}},{W_{0}^{c} = {\frac{\lambda}{n + \lambda} + \left( {1 - a^{2} + \gamma} \right)}},{W_{j}^{m} = {W_{j}^{c} = \frac{\lambda}{2\left( {n + \lambda} \right)}}},{j = 1},2,\ldots \;,{2n}$

where λ=(α²−1)n is a scaling parameter, n is the dimension of the statex_(α)(k), α is a parameter that determines the spread of the sigmapoints around x(k−1|k−1), which is generally assigned a small positivevalue, 10⁻³ for example, γ is a parameter used to incorporate a prioriknowledge on the distribution of x: for a Gaussian distribution, γ=2 isoptimal.

At the time k−1, the following selection of sigma points (set of pointsencoding exactly the mean and covariance information s considered):

${{x_{0}\left( {k - 1} \right)} = {{\hat{x}}_{a}\left( {k - 1} \middle| {k - 1} \right)}},{{x_{i}\left( {k - 1} \right)} = {{{\hat{x}}_{a}\left( {k - 1} \middle| {k - 1} \right)} + {\sqrt{{n + \lambda}\;}{S_{i}\left( {k - 1} \right)}}}},{i = 1},2,\ldots \mspace{14mu},n$${{x_{i + n}\left( {k - 1} \right)} = {{{\hat{x}}_{a}\left( {k - 1} \middle| {k - 1} \right)} - {\sqrt{{n + \lambda}\;}{S_{i}\left( {k - 1} \right)}}}},{i = 1},2,\ldots \mspace{14mu},n$

where S_(i)(k−1) is the i-th column of the matrix square root ofP_(x)(k−1|k−1), that is

P_(x)(k|k − 1) = S(k − 1)^(T)S(k − 1).

2. Prediction update

Each sigma point is propagated through the non-linear model representingthe evolution of the states:

x̂_(j)(k|k − 1) = A_(a)x_(j)(k − 1) + f_(a)(x_(j)(k − 1)) + B_(a)u(k − 1), j = 0, 1, …  , 2n

The mean and the covariance of {circumflex over (x)}_(α)(k|k−1), theprediction od x_(α)(k|K−1) are calculated as

$\mspace{20mu} {{{\hat{x}\left( k \middle| {k - 1} \right)} = {\sum\limits_{j = 0}^{2n}{W_{j}^{m}{{\hat{x}}_{j}\left( k \middle| {k - 1} \right)}}}},{{P_{x}\left( k \middle| {k - 1} \right)} = {{\sum\limits_{j = 0}^{2n}{{W_{j}^{\sigma}\left( {{{\hat{x}}_{j}\left( k \middle| {k - 1} \right)} - {x_{a}\left( k \middle| {k - 1} \right)}} \right)}\left( {{{\hat{x}}_{j}\left( k \middle| {k - 1} \right)} - {x_{a}\left( k \middle| {k - 1} \right)}} \right)^{T}}} + Q}}}$

The predicted states {circumflex over (x)}_(j)(k|k−1) are used in theoutput state equation, which yields:

y_(j)(k|k − 1) = C_(a)x_(j)(k|−k − 1)

The mean and the covariance of ŷ(k|k−1) are calculated as

$\mspace{20mu} {{{\hat{y}\left( k \middle| {k - 1} \right)} = {\sum_{j = 0}^{2n}{W_{j}^{m}{\hat{y}\left( k \middle| {k - 1} \right)}}}},{{P_{y}\left( k \middle| {k - 1} \right)} = {{\sum_{j = 0}^{2n}{{W_{j}^{c}\left( {{{\hat{y}}_{j}\left( k \middle| {k - 1} \right)} - {y\left( k \middle| {k - 1} \right)}} \right)}\left( {{\hat{y}}_{j}\left( k \middle| {k - 1} \right)} \right)^{T}}} + R}}}$

while the cross-covariance between {circumflex over (x)}_(α)(k|k−1) andŷ(k|k−1) is:

${P_{xy}\left( {k{k - 1}} \right)} = {\sum\limits_{j = 0}^{2n}\; {{W_{j}^{c}\left( {{{\hat{x}}_{j}\left( {k{k - 1}} \right)} - {x_{a}\left( {k{k - 1}} \right)}} \right)}\left( {{{\hat{y}}_{j}\left( {k{k - 1}} \right)} - {y\left( {k{k - 1}} \right)}} \right)^{T}}}$

3. Update from the measurements

As in the linear Kalman filter, the final state estimation is obtainedby correcting the prediction with a feedback on the output predictionerror (measured):

x̂(k) = x̂_(a)(k|k − 1) + K(ŷ(k) − ŷ(k|k − 1))

where gain K is given by:

K = P_(xy)(k|k − 1)P_(y)(k|k − 1)⁻¹

The a posterior covariance of the estimation is updated With the formulaas follows:

P_(x)(k|k) = P_(x)(k|k − 1) − KP_(y)(k|k − 1)K^(T)

The following rules are used to select matrices P₀ (covariance of theinitial state) and R:

-   -   if the initial state {circumflex over (x)}_(α)(k) at the time        k=0 is known, that is m(0)        x_(α)(0), then P₀ ⁻¹ is large.    -   if there is much noise in the measurements γ(k), then R is        small.

It is much more complex to select Q. It is generally chosen in diagonalform, as follows:

$Q = \begin{bmatrix}Q_{x} & 0 \\0 & Q_{m}\end{bmatrix}$

with Q_(m)>>Q_(x).

The method based on the UKF and the random walk model of the wave forceaccording to this implementation can be summarized as follows:

-   -   initializing k=0, state vector {circumflex over        (x)}_(α)(0|0)=m(0) and the state of the covariance matrix        P(0|0)=        .    -   at any time k:        -   using:        -   the measurements of the mobile part's position and velocity            γ(k)=[δ(k)            (k)] and of the force exerted by the PTO on the mobile part            u(k)        -   the results of the estimations of the previous step            {circumflex over (x)}_(α)(k−1|K−1), P(k−1|k−1)        -   parameters Q, R (covariance matrices)        -   determining the excitation force τ_(w) exerted by the            incident waves on the mobile part, denoted by            _(w)(k) for this embodiment, by carrying out the following            steps;        -   applying the three steps of the unscented Kalman fitter            algorithm to obtain {circumflex over (x)}_(α)(k|K), P(k|k)            as described above, the complete state being thus estimated            with its covariance matrix:

x̂(k) = x̂_(a)(k|k − 1) + K(ŷ(k) − ŷ(k|k − 1))K = P_(xy)(k|k − 1)P_(y)(k|k − 1)⁻¹P_(x)(k|k) = P_(x)(k|k − 1) − KP_(y)(k|k − 1)K^(T)

-   -     with x_(α)(k) being the unknown state vector

${{x_{a}(k)} = \begin{bmatrix}{x(k)} \\{\tau_{w}(k)}\end{bmatrix}},$

-   -             matrices of the state-space realization, P being the covariance        matrix of the state vector, Q and R being calibration matrices,        and        -   extracting the last component of the estimated state vector,            that is the excitation force τ_(w) exerted by the incident            waves on the mobile part, by an equation of the form:            ŵ(k)=[0 1]{circumflex over (x)}_(α)(k|K).

7—Control of the Wave Energy System

This is an optional step. In this step, the wave energy system iscontrolled to account for the excitation force exerted by the incidentwaves. It is thus possible to drive the wave energy system in order tooptimize the recovered energy.

This step can control the mobile part of the wave energy system, usingfor example an electrical, pneumatic or hydraulical conversion machinereferred to as PTO (Power Take-Off) system. This PTO system influencesthe motion of the mobile part and allows the mechanical energy to betransferred to the electrical, pneumatic or hydraulical network. Modelpredictive control (MPC) is an example of a method of controlling waveenergy systems.

Comparative Examples

The features and advantages of the method according to the inventionwill be clear from reading the comparative examples hereafter.

For these examples, a flap wave energy system is used as illustrated inFIG. 1, at small scale for basin testing.

The hydrodynamic parameters of the wave energy system are described inTable 1.

TABLE 1 Hydrodynamic parameters of the wave energy system Hydrodynamicparamenters Total moment of inertia I_(eq) + I

0.2753 kg · m² Hydrostatic stiffness K 118.25N · m · rad⁻¹ coefficientDrag coefficient β 2N · m · rad⁻¹ · s²

indicates data missing or illegible when filed

Furthermore, for the state-space representation of the wave energysystem, the following matrices are used:

${A_{r} = \begin{bmatrix}{{- 5}\text{,}1899} & {{- 36}\text{,}7434} \\1 & 0\end{bmatrix}},{B_{r} = \begin{bmatrix}1 \\0\end{bmatrix}},{C_{r} = \left\lbrack {13\text{,}2612\mspace{14mu} 0} \right\rbrack},{D_{r} = 0}$

First, the excitation force exerted by the waves on the mobile part isdetermined for this wave energy system by use of the method according tothe prior art described in patent application FR-3,049,989, using alinear Kalman filter after linearization of the drag term. FIG. 2illustrates the curves of the excitation torque τ_(w) in Nm exerted bythe waves on the mobile part as a function of time T in s, for a regularwave. Curve REF in solid line corresponds to the reference torque reallyexerted and curve EST in dotted line corresponds to the torque estimatedwith this method of the prior art. It is observed that the methodaccording to the prior art with the linear Kalman filter gives impreciseresults, with an imprecision that increases with the wave amplitude.

Second, the excitation force exerted by the waves on the mobile part isdetermined for this wave energy system by use of the method according tothe invention (using an unscented Kalman filter). FIG. 3 illustrates thecurves of the excitation torque τ_(w) in Nm exerted by the waves on themobile part as a function of time T in s, for a regular wave. Curve REFin solid line corresponds to the reference torque really exerted andcurve EST in dotted line corresponds to the torque estimated with themethod according to the invention. Curves REF and EST are superposed,which shows that the method according to the invention allows precisedetermination of the excitation torque exerted by the waves on themobile part.

Thirdly, the second test of FIG. 2 is repeated for irregular waves. FIG.4 illustrates the curves of the excitation torque τ_(w) in Nm exerted bythe waves on the mobile part as a function of time T in s, for anirregular wave. Curve REF in solid line corresponds to the referencetorque really exerted and curve EST in dotted line corresponds to thetorque estimated with the method according to the invention. Curves REFand EST are superposed, which shows that the method according to theinvention allows precise determination of the excitation torque exertedby the waves on the mobile part, even for irregular waves.

The method according to the invention is therefore well suited forprecisely determining the excitation force exerted by the waves on themobile part, in particular of a non-linear wave energy system.

1-10. (canceled)
 11. A method of determining an excitation force exerted by incident waves on a mobile part of a wave energy system, the wave energy system converting the wave energy to electrical energy through the mobile part cooperating with a conversion machine, the mobile part moving with respect to the conversion machine under action of the waves, comprising: a) measuring position and velocity of the mobile part; b) measuring force exerted by the conversion machine on the mobile part; c) constructing a model of radiation force exerted on the mobile part, the radiation force model relating the radiation force to the velocity of the mobile part; d) constructing a model of drag force exerted on the mobile part, the drag force model relating the drag force to the velocity of the mobile part; e) constructing a dynamic model of the wave energy system relating the excitation force exerted by the incident waves on the mobile part to position of the mobile part, to velocity of the mobile part, to force exerted by the conversion machine on the mobile part, to radiation force exerted on the mobile part and to the drag force exerted on the mobile part; and f) determining the excitation force exerted by the incident waves on the mobile part using the dynamic model, the radiation force model, the drag force model, the measured position and the measured velocity, and the measured force exerted by the conversion machine on the mobile part, and using an unscented Kalman filter constructed from a random walk model of the excitation force exerted by the incident waves on the mobile part.
 12. A method as claimed in claim 11, wherein the dynamic model of the wave energy system is constructed by use of an equation:

(t)=τ_(hd)(t)+τ_(rad)(t)+τ_(d)(t)+τ_(w)(t)−u(t), with

being a total moment of inertia of the mobile part, δ(t) being an angle of rotation of the mobile part with respect to an equilibrium position, with

(t) being an angular acceleration of the mobile part and

(t) being an angular velocity of the mobile part, τ_(hd)(t) being a hydrostatic restoring moment, τ_(rad)(t) being a radiation moment, τ_(d)(t) being a drag moment, τ_(w)(t) being a wave excitation moment and u(t) being a moment exerted by the conversion machine on the mobile part.
 13. A method as claimed in claim 12, wherein hydrostatic restoring moment τ_(hd)(t) is determined by use of a formula: τ_(hd)(t) = −K δ(t) where K is a hydrostatic stiffness coefficient.
 14. A method as claimed in claim 11, wherein the radiation force model is constructed by use of an equation: τ_(rad)(t)=−

(t)−τ_(r)(t), with

being an added moment of inertia at an infinitely high frequency, and τ_(r)(t)=

h(t−s)

(s)ds=h(t)

(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and

(t) being an angular velocity of the mobile part.
 15. A method as claimed in claim 12, wherein the radiation force model is constructed by use of an equation: τ_(rad)(t)=−

(t)−τ_(r)(t), with

being an added moment of inertia at an infinitely high frequency, and τ_(r)(t)=

h(t−s)

(s)ds=h(t)

(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and

(t) being an angular velocity of the mobile part.
 16. A method as claimed in claim 13, wherein the radiation force model is constructed by use of an equation: τ_(rad)(t)=−

(t)−τ_(r)(t), with

being an added moment of inertia at an infinitely high frequency, and τ_(r)(t)=

h(t−s)

(s)ds=h(t)

(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and

(t) being an angular velocity of the mobile part.
 17. A method as claimed in claim 14, wherein the radiation force model is constructed by use of an equation: τ_(rad)(t)=−

(t)−τ_(r)(t), with

being an added moment of inertia at an infinitely high frequency, and τ_(r)(t)=

h(t−s)

(s)ds=h(t)

(t), with h being an impulse response relating the velocity of the mobile part to a radiation damping and

(t) being an angular velocity of the mobile part.
 18. A method as claimed in claim 11, wherein the drag force model is constructed by use of an equation: ${\tau_{d}(t)} = {\beta {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$ where β is the drag coefficient and

(t) is the angular velocity of the mobile part.
 19. A method as claimed in claim 12, wherein the drag force model is constructed by use of an equation: ${\tau_{d}(t)} = {\beta {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$ where β is the drag coefficient and

(t) is the angular velocity of the mobile part.
 20. A method as claimed in claim 13, wherein the drag force model is constructed by use of an equation: ${\tau_{d}(t)} = {\beta {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$ where β is the drag coefficient and

(t) is the angular velocity of the mobile part.
 21. A method as claimed in claim 14, wherein the drag force model is constructed by use of an equation: ${\tau_{d}(t)} = {\beta {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$ where β is the drag coefficient and

(t) is the angular velocity of the mobile part.
 22. A method as claimed in claim 15, wherein the drag force model is constructed by use of an equation: ${\tau_{d}(t)} = {\beta {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$ where β is the drag coefficient and

(t) is the angular velocity of the mobile part.
 23. A method as claimed in claim 16, wherein the drag force model is constructed by use of an equation: ${\tau_{d}(t)} = {\beta {\overset{.}{\delta}(t)}{{\overset{.}{\delta}(t)}}}$ where β is the drag coefficient and

(t) is the angular velocity of the mobile part.
 24. A method as claimed claim 11, wherein the drag force model F_(d) is expressed by an equation: ${F_{d}(t)} = {\beta {\overset{.}{z}(t)}{{\overset{.}{z}(t)}}}$ where β is the drag coefficient and

being velocity of the mobile part.
 25. A method as claimed in claim 11, wherein the excitation force τ_(w) exerted by the incident waves on the mobile part is determined by carrying out steps of: i) initializing k=0, state vector {circumflex over (x)}_(α)(0|0)=m(0) and a state of the covariance matrix P(0|0)=P₀, ii) at any time k, acquiring the position and velocity measurements of the mobile part γ(k)=[δ(k)

(k)] and the measurement of the force exerted by the conversion machine on the mobile part u(k), and iii) at any time k, determining the excitation force exerted by the incident waves on the mobile part

_(w)(k) by use of equations as follows:      τ̂_(w)(k) = [0  1]x̂_(a)(kk)      P_(x)(kk) = P_(x)(kk − 1) − KP_(y)(kk − 1)K^(T)      K = P_(xy)(kk − 1)P_(y)(kk − 1)⁻¹      x̂(k) = x̂_(a)(kk − 1) + K(ŷ(k) − ŷ(kk − 1)) ${P_{xy}\left( {k{k - 1}} \right)} = {\sum\limits_{j = 0}^{2n}\; {{W_{j}^{c}\left( {{{\hat{x}}_{j}\left( {k{k - 1}} \right)} - {x_{a}\left( {k{k - 1}} \right)}} \right)}\left( {{{\hat{y}}_{j}\left( {k{k - 1}} \right)} - {y\left( {k{k - 1}} \right)}} \right)^{T}}}$ $\mspace{76mu} {{{\hat{y}\left( {k{k - 1}} \right)} = {\sum\limits_{j = 0}^{2n}\; {W_{j}^{m}{{\hat{y}}_{j}\left( {k{k - 1}} \right)}}}},{{P_{y}\left( {k{k - 1}} \right)} = {{\sum\limits_{j = 0}^{2n}\; {{W_{j}^{c}\left( {{{\hat{y}}_{j}\left( {k{k - 1}} \right)} - {y\left( {k{k - 1}} \right)}} \right)}\left( {{{\hat{y}}_{j}\left( {k{k - 1}} \right)} - {y\left( {k{k - 1}} \right)}} \right)^{T}}} + R}}}$      ŷ_(j)(kk − 1) = C_(a)x̂_(j)(kk − 1) $\mspace{76mu} {{{\hat{x}\left( {k{k - 1}} \right)} = {\sum\limits_{j = 0}^{2n}\; {W_{j}^{m}{{\hat{x}}_{j}\left( {k{k - 1}} \right)}}}},{{P_{x}\left( {k{k - 1}} \right)} = {{\sum\limits_{j = 0}^{2n}\; {{W_{j}^{c}\left( {{{\hat{x}}_{j}\left( {k{k - 1}} \right)} - {x_{a}\left( {k{k - 1}} \right)}} \right)}\left( {{{\hat{x}}_{j}\left( {k{k - 1}} \right)} - {x_{a}\left( {k{k - 1}} \right)}} \right)^{T}}} + Q}}}$ x̂_(j)(kk − 1) = A_(a)x_(j)(k − 1) + f_(a)(x_(j)(k − 1)) + B_(a)u(k − 1), j = 0, 1, … , 2n $\mspace{76mu} {{{x_{0}\left( {k - 1} \right)} = {{\hat{x}}_{a}\left( {{k - 1}{k - 1}} \right)}},\mspace{76mu} {{x_{i}\left( {k - 1} \right)} = {{{\hat{x}}_{a}\left( {{k - 1}{k - 1}} \right)} + {\sqrt{n + \lambda}{S_{i}\left( {k - 1} \right)}}}},{i = 1},2,\ldots \;,n}$ $\mspace{76mu} {{{x_{i + n}\left( {k - 1} \right)} = {{{\hat{x}}_{a}\left( {{k - 1}{k - 1}} \right)} - {\sqrt{n + \lambda}{S_{i}\left( {k - 1} \right)}}}},{i = 1},2,\ldots \;,n}$ where S_(i)(k−1) is an i-th column of a matrix square root of P_(x)(k−1|K−1), is P_(x)(k|k − 1) = S(k − 1)^(T)S(k − 1) x_(α)(k) is an unknown state vector ${{x_{a}(k)} = \begin{bmatrix} {x(k)} \\ {\tau_{w}(k)} \end{bmatrix}},$

and matrices of the state-space realization, P is a covariance matrix of an unknown state vector, and Q and R are calibration matrices.
 26. A method as claimed in claim 11, wherein the wave energy system is controlled depending on the determined excitation force τ_(w) exerted by the incident waves on the mobile part.
 27. A method as claimed in claim 11, wherein the wave energy system is one of a wave energy system including submerged pressure differential converter device one of a surface floating rotating mass device or a flap device wave energy system.
 28. A wave energy system, comprising a control, including a computer control which implements a method as claimed in claim
 11. 